Pour ce TP, nous allons utiliser les librairies de R suivantes :

library(reticulate)
library(ggplot2)
library(gridExtra)
library(plotly)
library(FactoMineR)
library(factoextra)
library(dplyr)
library(klaR)
library(reshape2)
library(mclust)
library(Rmixmod)
library(cluster)

et les librairies de Python suivantes :

import numpy as np
import pandas as pd
import plotly.express as px
import matplotlib.pyplot as plt

import prince
from kmodes.kmodes import KModes

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

from yellowbrick.cluster import SilhouetteVisualizer

import gower
from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.cluster import AgglomerativeClustering

from sklearn.metrics.cluster import adjusted_rand_score

L’objectif de ce TP est d’étudier quelques stratégies pour faire du clustering d’individus décrits par des variables qualitatives.

1 Jeu de données et statistiques descriptives

1.1 Description des données mushroom

Nous allons nous intéresser au jeu de données de champignons disponible sur le site de l’UCI : mushroom. Pour limiter le temps de calcul et faciliter l’interprétation durant le TP, nous allons considérer que \(800\) individus (\(400\) de chaque classe) parmi les \(8124\) disponibles et l’on se restreint aux variables suivantes :

    1. Class: edible (e) or poisonous (p)
    1. odor: almond=a,anise=l,creosote=c,fishy=y,foul=f, musty=m,none=n,pungent=p,spicy=s
    1. gill-color: black=k,brown=n,buff=b,chocolate=h,gray=g, green=r,orange=o,pink=p,purple=u,red=e, white=w,yellow=y
    1. stalk-root: bulbous=b,club=c,cup=u,equal=e, rhizomorphs=z,rooted=r,missing=?
    1. stalk-color-above-ring: brown=n,buff=b,cinnamon=c,gray=g,orange=o, pink=p,red=e,white=w,yellow=y
    1. stalk-color-below-ring: brown=n,buff=b,cinnamon=c,gray=g,orange=o, pink=p,red=e,white=w,yellow=y
    1. ring-number: none=n,one=o,two=t
    1. ring-type: cobwebby=c,evanescent=e,flaring=f,large=l, none=n,pendant=p,sheathing=s,zone=z
    1. spore-print-color: black=k,brown=n,buff=b,chocolate=h,green=r, orange=o,purple=u,white=w,yellow=y
    1. population: abundant=a,clustered=c,numerous=n, scattered=s,several=v,solitary=y
    1. habitat: grasses=g,leaves=l,meadows=m,paths=p, urban=u,waste=w,woods=d

Dans la suite, on note la matrice des données \[\mathbf{X}=(x_{ij})_{\begin{array}{l} i=1,\ldots,n\\ j=1,\ldots,p \end{array}}\] où \(x_{ij}\) est la modalité prise par le \(i\)-ème individu pour la \(j\)-ème variable qualitative.

Question : Chargez le jeu de données disponible dans le fichier mushroomTP.csv déposé sur la page moodle du cours. Conservez à part la variable Class.

Data <- read.table("...")
Classep <- ...
Data <- ....  # jeu de données sans la variable Class

Correction :

Data <- read.table("mushroomTP.csv", header = T, sep = ",", stringsAsFactors = T)
classesp <- Data[, 1]
Data <- Data[, -1]

Afin de plus facilement repérer les variables et leur modalité, on modifie le jeu de données :

for (j in 1:ncol(Data)) {
    if (j < 10) {
        Data[, j] <- as.factor(paste("0", j, "-", Data[, j], sep = ""))
    } else {
        Data[, j] <- as.factor(paste(j, "-", Data[, j], sep = ""))
    }
}

Question : Vérifiez la taille du jeu de données et faites un summary pour contrôler le contenu du jeu de données Data.

# A FAIRE

Correction :

dim(Data)
[1] 800  10
summary(Data)
      odor       gill.color  stalk.root stalk.color.above.ring
 01-n   :333   02-b   :188   03-?:260   04-w   :424           
 01-f   :204   02-p   :148   03-b:359   04-p   :198           
 01-y   : 69   02-w   :113   03-c: 54   04-g   : 58           
 01-s   : 63   02-n   : 91   03-e:105   04-b   : 47           
 01-a   : 42   02-g   : 73   03-r: 22   04-n   : 39           
 01-l   : 38   02-h   : 73              04-o   : 17           
 (Other): 51   (Other):114              (Other): 17           
 stalk.color.below.ring ring.number ring.type  spore.print.color population
 05-w   :434            06-n:  6    07-e:278   08-w   :254       09-a: 30  
 05-p   :194            06-o:734    07-f:  6   08-k   :191       09-c: 35  
 05-g   : 48            06-t: 60    07-l:127   08-n   :178       09-n: 31  
 05-b   : 47                        07-n:  6   08-h   :154       09-s:120  
 05-n   : 44                        07-p:383   08-r   :  9       09-v:413  
 05-o   : 17                                   08-o   :  5       09-y:171  
 (Other): 16                                   (Other):  9                 
 habitat   
 10-d:300  
 10-g:196  
 10-l: 93  
 10-m: 30  
 10-p:129  
 10-u: 33  
 10-w: 19

De même on lit le jeu de données sous Python :

data_full = pd.read_csv('mushroomTP.csv')
target = data_full[['class']]
Datapy = data_full.drop(['class'],axis=1)

1.2 Analyse des correspondances multiples

Question : Faites une analyse des correspondances multiples avec R. Vous conserverez les coordonnées des projections des individus pour la suite du TP.

library(FactoMineR)
library(factoextra)
resacm <- MCA(....)
fviz_screeplot(resacm, choice = "eigenvalue")
fviz_mca(...)

Correction :

library(FactoMineR)
resacm <- MCA(Data, ncp = 10, graph = FALSE)
fviz_screeplot(resacm, choice = "eigenvalue")

ggplotly(fviz_mca(resacm, geom = "point", geom.var = "text"))
ggplotly(fviz_mca(resacm, geom = "point", geom.var = "text", axes = c(1, 3)))
fviz_mca_ind(resacm, geom = "point", habillage = classesp)

fviz_mca_var(resacm, col.var = "contrib") + scale_color_gradient2(low = "white",
    mid = "blue", high = "red", midpoint = 2) + theme_minimal()

fviz_contrib(resacm, choice = "var", axes = 1, top = 10)

Question : Faites également une analyse des correspondances multiples avec Python.

# A COMPLETER
import prince 
mca=prince.MCA(....)
mca = mca.fit(Datapy)

Correction :

import prince
mca = prince.MCA(
     n_components=10,
     n_iter=3,
     copy=True,
     check_input=True,
     engine='auto',
     random_state=42
 ).fit(Datapy)

mca_df = pd.DataFrame({
    "Dim1" : mca.row_coordinates(Datapy)[0], 
    "Dim2" : mca.row_coordinates(Datapy)[1],
    "Class": data_full["class"],
    "Var1" : mca.column_coordinates(Datapy)[0],
    "Var2" : mca.column_coordinates(Datapy)[1]
})

fig=px.scatter(mca_df,x="Dim1",y="Dim2",color="Class")
fig.show()

1.3 Ingrédients pour la suite

On se munit également d’une fonction auxiliaire pour visualiser les résultats du clustering pour la suite. Cette fonction prend en argument un clustering clust, le jeu de données Data et un vecteur J contenant les indices des variables d’intérêt. Il renvoie la fréquence des modalités des variables choisies vis-à-vis des classes de clust.

# J indice des noms de variables

heatm <- function(clust, Data, J) {
    library(dplyr)
    Dataaux <- data.frame(id.s = c(1:nrow(Data)), Data)
    aux <- cbind(Dataaux, clust)
    aux.long <- melt(data.frame(lapply(aux, as.character)), stringsAsFactors = FALSE,
        id = c("id.s", "clust"), factorsAsStrings = T)
    # Effectifs
    aux.long.q <- aux.long %>%
        group_by(clust, variable, value) %>%
        mutate(count = n_distinct(id.s)) %>%
        distinct(clust, variable, value, count)
    # avec fréquences
    aux.long.p <- aux.long.q %>%
        group_by(clust, variable) %>%
        mutate(perc = count/sum(count)) %>%
        arrange(clust)

    Lev <- NULL
    for (j in 1:ncol(Data)) Lev <- c(Lev, levels(Data[, j]))

    Jaux <- NULL
    for (j in 1:length(J)) {
        Jaux <- c(Jaux, which(aux.long.p$variable == colnames(Data)[J[j]]))
    }

    gaux <- ggplot(aux.long.p[Jaux, ], aes(x = clust, y = value)) + geom_tile(aes(fill = perc)) +
        scale_fill_gradient2(low = "white", mid = "blue", high = "red") + theme_minimal()

    return(list(gaux = gaux, eff = aux.long.q, freq = aux.long.p))
}

Mode <- function(x) {
    ux <- unique(x)
    ux[which.max(tabulate(match(x, ux)))]
}

barplotClus <- function(clust, Data, J) {
    aux.long.p <- heatm(clust, Data, J)$freq
    # ux<-unique(aux.long.p$variable)
    for (j in J) {
        p <- ggplot(aux.long.p[which(aux.long.p$variable == colnames(Data)[j]), ],
            aes(x = clust, y = perc, fill = value)) + geom_bar(stat = "identity")
        print(p)
    }

}

2 Méthode des Kmodes

Dans cette partie, nous allons utiliser la méthode des Kmodes introduite par Huang (1998). Rappelons que cette méthode est une extension des Kmeans dans le cas des données qualitatives. Les modifications par rapport aux Kmeans sont

  • le changement de distance : on utilise la dissimilarité basée sur l’appariement simple

\[ d(\mathbf{x}_i,\mathbf{x}_\ell) = \underset{j=1}{\stackrel{p}{\sum}}\ \mathbb{1}_{x_{ij}\neq x_{\ell j}} \]

  • le centre d’une classe est calculé en fonction des fréquences des modalités majoritaires présentes dans cette classe: pour la classe \(\mathcal{C}_k\),

\[ \mathbf{m}_k=(m_{k1},\ldots,m_{kp}) \textrm{ avec } m_{kj}= \underset{u_1,\ldots,u_{s_j}}{\mbox{argmax}}\ \underset{i\in\mathcal C_k}{\sum}\ \mathbb{1}_{x_{ij}= u_{s_j}} \]

Question : A l’aide de la fonction kmodes()de la librairie klaR, déterminez une classification en \(K=6\) classes des champignons. Comparez avec la classe de l’espèce. Visualisez la classification obtenue. Vous pouvez vous aider des fonctions auxiliaires pour interpréter la classification.

library(klaR)
reskmodes <- kmodes(....)

Correction :

library(klaR)
K <- 6
clkmodes <- kmodes(Data, K, iter.max = 100, weight = FALSE)
table(classesp, clkmodes$cluster)
        
classesp   1   2   3   4   5   6
       e 134 101  29   0 136   0
       p   8  38 194  60  32  68
clkmodes$modes
  odor gill.color stalk.root stalk.color.above.ring stalk.color.below.ring
1 01-n       02-n       03-b                   04-w                   05-w
2 01-n       02-p       03-b                   04-w                   05-w
3 01-y       02-b       03-?                   04-w                   05-p
4 01-f       02-h       03-b                   04-p                   05-b
5 01-n       02-p       03-e                   04-w                   05-w
6 01-f       02-g       03-b                   04-b                   05-n
  ring.number ring.type spore.print.color population habitat
1        06-o      07-p              08-n       09-y    10-d
2        06-o      07-p              08-k       09-v    10-d
3        06-o      07-e              08-w       09-v    10-l
4        06-o      07-l              08-h       09-v    10-d
5        06-o      07-p              08-k       09-s    10-g
6        06-o      07-l              08-h       09-v    10-p
fviz_mca_ind(resacm, habillage = as.factor(clkmodes$cluster), geom = "point") + labs(title = "Kmodes")

fviz_mca_ind(resacm, habillage = as.factor(clkmodes$cluster), geom = "point", axes = c(1,
    3)) + labs(title = "Kmodes")

fviz_mca(resacm, habillage = as.factor(clkmodes$cluster), geom.ind = c("point"),
    geom.var = c("point", "text"))

# heatm(clust=clkmodes$cluster,Data=Data,J=c(1:10)$gaux
barplotClus(clkmodes$cluster, Data, J = c(1:10))

Question : Pour déterminer le nombre de classes, la méthode du coude peut être utilisée en remplaçant l’inertie intra-classe par le critère “Within Cluster Difference”

\[ WCD(K) = \underset{k=1}{\stackrel{K}{\sum}} \underset{i\in\mathcal C_k}{\sum}\ d(x_i,m_k) \]

où \(d(.,.)\) est l’appariement simple et \(m_k\) est le centre de la classe \(\mathcal C_k\).

Tracez la courbe \(K\mapsto WCD(K)\) pour déterminer le nombre de classes optimal. Vous pouvez vous aider des sorties de la fonction kmodes().

# A COMPLETER
WithinDiff <- NULL
Kmax <- 10
Clust <- matrix(0, nrow = nrow(Data), ncol = Kmax)
for (k in 1:Kmax) {
    aux <- kmodes(.........)
    WithinDiff <- c(WithinDiff, ..........)
    Clust[, k] <- aux$cluster
}

auxdf <- data.frame(NbCluster = 1:Kmax, WithinDiff = WithinDiff)
ggplot(auxdf, aes(x = NbCluster, y = WithinDiff)) + geom_point() + geom_line()

Correction :

WithinDiff <- NULL
Kmax <- 10
Clust <- matrix(0, nrow = nrow(Data), ncol = Kmax)
for (k in 1:Kmax) {
    aux <- kmodes(Data, k, iter.max = 100, weight = FALSE)
    WithinDiff <- c(WithinDiff, sum(aux$withindiff))
    Clust[, k] <- aux$cluster
}

auxdf <- data.frame(NbCluster = 1:Kmax, WithinDiff = WithinDiff)
ggplot(auxdf, aes(x = NbCluster, y = WithinDiff)) + geom_point() + geom_line()

Question : Reprenez les questions précédentes avec Python. Vous pouvez utiliser la fonction KModes.

from kmodes.kmodes import KModes

cost = []
for clustk in range(1, 10):
  kmodes = KModes(n_jobs = -1, n_clusters = clustk, init ='Huang',random_state = 0)
  kmodes.fit_predict(Datapy);
  cost.append(kmodes.cost_);

# A COMPLETER

Correction :

import numpy as np
from kmodes.kmodes import KModes

# Cout
cost = []
for clustk in range(1, 10):
  kmodes = KModes(n_jobs = -1, n_clusters = clustk, init ='Huang',random_state = 0)
  res=kmodes.fit_predict(Datapy);
  cost.append(kmodes.cost_);
        
# Converting the results into a dataframe and plotting them
df_cost = pd.DataFrame({'Cluster': range(1, 10), 'Cost': cost})

# Data viz
import plotly.express as px
fig = px.line(df_cost, x='Cluster', y='Cost',markers=True)
fig.show()
kmodes = KModes(n_jobs = -1, n_clusters = 6, init ='Huang',random_state = 0);
kmodes.fit(Datapy);
pd.DataFrame(kmodes.cluster_centroids_,columns=Datapy.columns)
  odor gill.color stalk.root  ... spore.print.color population habitat
0    n          n          b  ...                 k          v       d
1    n          w          ?  ...                 w          c       g
2    n          p          e  ...                 n          s       g
3    f          h          b  ...                 h          v       p
4    n          p          b  ...                 k          y       d
5    y          b          ?  ...                 w          v       l

[6 rows x 10 columns]
pd.crosstab(kmodes.labels_,"freq")
col_0  freq
row_0      
0       228
1        62
2       109
3       124
4        80
5       197

3 Méthodes des Kmeans sur coordonnées de l’ACM

Une seconde stratégie est de partir des coordonnées de l’analyse des correspondances multiples et d’utiliser un algorithme plus usuel sur données quantitatives. Dans cette section, on va appliquer l’algorithme des Kmeans.

Question : Appliquez l’algorithme des Kmeans sur les coordonnées de l’ACM. Pour la détermination du nombre de classes, vous pouvez utiliser l’évolution de l’inertie intra-classe et le critère silhouette.

# A COMPLETER

Correction :

set.seed(1234)
Kmax <- 15
kmeansclus <- matrix(0, nrow = nrow(Data), ncol = (Kmax - 1))
Iintra <- NULL
Silhou <- NULL
for (k in 2:Kmax) {
    resaux <- kmeans(resacm$ind$coord[, 1:5], centers = k, nstart = 10)
    kmeansclus[, (k - 1)] <- resaux$cluster
    Iintra <- c(Iintra, resaux$tot.withinss)
    aux <- silhouette(resaux$cluster, daisy(resacm$ind$coord))
    Silhou <- c(Silhou, mean(aux[, 3]))
}

df <- data.frame(K = 2:Kmax, Iintra = Iintra, Silhou = Silhou)
g1 <- ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point() + xlab("Nombre de classes") +
    ylab("Inertie intraclasse")
g2 <- ggplot(df, aes(x = K, y = Silhou)) + geom_line() + geom_point() + xlab("Nombre de classes") +
    ylab("Critère Silhouette")
grid.arrange(g1, g2, ncol = 2)

On retient par la suite \(8\) classes.

Question : Etudiez la classification retenue. Comparez avec la classification obtenue précédemment avec les Kmodes et la variable “Class” (champignon comestible ou toxique).

# A COMPLETER

Correction :

Kaux <- 8 - 1
table(kmeansclus[, Kaux])

  1   2   3   4   5   6   7   8 
 17 190 127  45   6 221 177  17 
aux <- silhouette(kmeansclus[, Kaux], daisy(resacm$ind$coord))
fviz_silhouette(aux) + theme(plot.title = element_text(size = 9))
  cluster size ave.sil.width
1       1   17          0.72
2       2  190          0.76
3       3  127          0.69
4       4   45          0.28
5       5    6          0.92
6       6  221          0.53
7       7  177          0.15
8       8   17          0.57

fviz_mca(resacm, geom.ind = "point", geom.var = c("point", "text"), habillage = as.factor(kmeansclus[,
    Kaux]))

fviz_mca_ind(resacm, geom = "point", habillage = as.factor(kmeansclus[, Kaux]), axes = c(1,
    2))

fviz_mca_ind(resacm, geom = "point", habillage = as.factor(kmeansclus[, Kaux]), axes = c(1,
    3))

# Comparaison avec Kmodes
table(clkmodes$cluster, kmeansclus[, Kaux])
   
      1   2   3   4   5   6   7   8
  1   0   0   0   1   0  94  37  10
  2   0   0   0  13   3 115   7   1
  3  17 190   0   4   2   5   0   5
  4   0   0  59   0   1   0   0   0
  5   0   0   0  27   0   7 133   1
  6   0   0  68   0   0   0   0   0
adjustedRandIndex(clkmodes$cluster, kmeansclus[, Kaux])
[1] 0.5784484
# Comparaison avec comestible/toxique
table(classesp, kmeansclus[, Kaux])
        
classesp   1   2   3   4   5   6   7   8
       e  17   0   0  34   0 186 146  17
       p   0 190 127  11   6  35  31   0
adjustedRandIndex(classesp, kmeansclus[, Kaux])
[1] 0.2779427
# heatm(clust=kmeansclus[,9],Data=Data,J=c(1:5))
barplotClus(kmeansclus[, Kaux], Data, J = c(1:10))

Question : Reprenez les questions précédentes avec Python.

# A completer

Correction :

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

mca = prince.MCA(n_components=5).fit(Datapy)

silhouette_coefficients = []
InertieIntra=[]
for k in range(3, 16):
  kmeans = KMeans(n_clusters=k,init = 'random', max_iter = 300, n_init = 10, random_state = 10);
  res=kmeans.fit(mca.row_coordinates(Datapy));
  score = silhouette_score(mca.row_coordinates(Datapy), kmeans.labels_);
  silhouette_coefficients.append(score);
  InertieIntra.append(kmeans.inertia_);

aux = pd.DataFrame({"K" : range(3, 16),
                    "Silhouette" : silhouette_coefficients,
                    "Intra" : InertieIntra})
fig = px.line(aux, x='K', y='Silhouette',labels={'x':'Number of clusters', 'y':'Silhouette'},markers=True)
fig.show() 
fig = px.line(aux, x='K', y='Intra',labels={'x':'Number of clusters', 'y':'Inertie Intra'},markers=True)
fig.show() 
from yellowbrick.cluster import SilhouetteVisualizer
kmeansopt = KMeans(6,init = 'k-means++', max_iter = 300, n_init = 10, random_state = 0);
visualizer = SilhouetteVisualizer(kmeansopt, colors='yellowbrick');
visualizer.fit(mca.row_coordinates(Datapy))     
SilhouetteVisualizer(ax=<AxesSubplot: >, colors='yellowbrick',
                     estimator=KMeans(n_clusters=6, random_state=0))
visualizer.show()        

kmeanslab=kmeansopt.predict(mca.row_coordinates(Datapy))
pd.crosstab(kmeanslab,"freq")
col_0  freq
row_0      
0       127
1       434
2         6
3        19
4       197
5        17
pd.crosstab(kmeanslab,kmodes.labels_)
col_0    0   1    2    3   4    5
row_0                            
0        0   0    0  124   3    0
1      215  33  109    0  77    0
2        0   6    0    0   0    0
3        0  19    0    0   0    0
4        0   0    0    0   0  197
5       13   4    0    0   0    0
adjusted_rand_score(kmeanslab,kmodes.labels_)
0.515926879468948
pd.crosstab(kmeanslab,data_full["class"])
class    e    p
row_0          
0        0  127
1      359   75
2        0    6
3       19    0
4        5  192
5       17    0
adjusted_rand_score(kmeanslab,data_full["class"])
0.4129657489670915

4 Méthode de classification hiérarchique

Question : Quelle classification par CAH pouvez-vous mettre en place pour traiter des données qualitatives ? Mettez en application votre proposition avec R et étudiez la classification retenue. Vous pourrez vous aider des fonctions daisy() de la librairie cluster, de la fonction hclust(), …

# A COMPLETER

Correction :

gower.dist <- daisy(Data, metric = c("gower"))
aggl <- hclust(gower.dist, method = "complete")
ggplot(data.frame(K = 1:20, Height = sort(aggl$height, decreasing = T)[1:20]), aes(x = K,
    y = Height)) + geom_line() + geom_point()

fviz_dend(aggl, show_labels = FALSE, k = 5)

clustaggl <- cutree(aggl, 5)
fviz_mca(resacm, geom.ind = "point", geom.var = c("point", "text"), habillage = as.factor(clustaggl))

table(clustaggl)
clustaggl
  1   2   3   4   5 
274 177 216 127   6 
table(clustaggl, classesp)
         classesp
clustaggl   e   p
        1 199  75
        2 177   0
        3  24 192
        4   0 127
        5   0   6
table(clustaggl, clkmodes$cluster)
         
clustaggl   1   2   3   4   5   6
        1  53  54   0   0 167   0
        2  89  82   5   0   1   0
        3   0   0 216   0   0   0
        4   0   0   0  59   0  68
        5   0   3   2   1   0   0
barplotClus(clust = clustaggl, Data, J = c(1:10))

Question : Reprenez la question avec python. Vous pouvez utiliser la librairie gower.

import gower
from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.cluster import AgglomerativeClustering

# A COMPLETER

Correction :

import gower
from scipy.spatial.distance import squareform
distance_matrix = gower.gower_matrix(Datapy);
condensed_dist_matrix = squareform(distance_matrix)

from scipy.cluster.hierarchy import dendrogram, linkage
Zd=linkage(condensed_dist_matrix,'complete');
res=dendrogram(Zd, orientation='top',distance_sort='descending',show_leaf_counts=False,no_labels=True);
plt.show()

from sklearn.cluster import AgglomerativeClustering
model = AgglomerativeClustering(n_clusters=9,affinity='precomputed',linkage='complete')
ClustK = model.fit_predict(distance_matrix)
pd.crosstab(ClustK,"freq")
col_0  freq
row_0      
0       148
1       196
2        20
3        41
4       162
5       127
6        83
7         6
8        17
pd.crosstab(ClustK,data_full["class"])
class    e    p
row_0          
0       82   66
1        5  191
2       19    1
3       32    9
4      162    0
5        0  127
6       83    0
7        0    6
8       17    0
pd.crosstab(ClustK,kmodes.labels_)
col_0    0   1   2    3   4    5
row_0                           
0      106   1  41    0   0    0
1        0   0   0    0   0  196
2        0  19   0    0   0    1
3        9  32   0    0   0    0
4       85   0   0    0  77    0
5        0   0   0  124   3    0
6       15   0  68    0   0    0
7        0   6   0    0   0    0
8       13   4   0    0   0    0
adjusted_rand_score(ClustK,kmodes.labels_)
0.6528431966802858
mca_df = pd.DataFrame({
    "Dim1" : mca.row_coordinates(Datapy)[0], 
    "Dim2" : mca.row_coordinates(Datapy)[1],
    "Dim3" : mca.row_coordinates(Datapy)[2],
    "Clust" : pd.Categorical(ClustK)
})

fig=px.scatter(mca_df,x="Dim1",y="Dim2",color="Clust")
fig.show()
fig=px.scatter(mca_df,x="Dim1",y="Dim3",color="Clust")
fig.show()

5 Méthode de classification par des mélanges

Pour cette dernière partie du TP, nous allons utiliser une stratégie de classification via des modèles de mélange. Rappelons que l’on a pour données \(\mathbf{X}=(x_1,\ldots,x_n)\) avec \(x_i\) décrit par \(p\) variables catégorielles, chacune avec \(m_j\) modalités. Dans ce cas de variables qualitatives, on peut considérer des distributions multinomiales par variables et par composantes. Pour écrire la distribution de mélange, on commence par transformer l’écriture des données de la façon suivante :

\[x_i=(x_{i1},\ldots,x_{ip})\rightsquigarrow (x_i^{jh}; j=1,\ldots,p, h=1,\ldots,m_j)\] avec \[ x_i^{\,jh}=\left\{\begin{array}{l l}1 & \textrm{ si } i \textrm{ prend la modalité } h \textrm{ pour la variable }j\\ 0 & \textrm{ sinon.}\end{array}\right. \] Les densités de mélange considérées sont de la forme \[f(.|\theta_K)=\underset{k=1}{\stackrel{K}{\sum}}\pi_k f_k(x_i|\boldsymbol{\alpha}_k)\] avec

  • \(f_k(x_i|\boldsymbol{\alpha}_k)=\underset{j=1}{\stackrel{p}{\prod}}\underset{h=1}{\stackrel{m_j}{\prod}}\left(\alpha_k^{\,jh}\right)^{x_i^{jh}}\)
  • \(\boldsymbol{\alpha}_k=(\alpha_k^{\,jh}; j=1,\ldots,p, h=1,\ldots,m_j)\) avec \(\underset{h=1}{\stackrel{m_j}{\sum}}\alpha_k^{jh}=1\)
  • \(\theta_k=(\pi_1,\ldots,\pi_K,\boldsymbol{\alpha}_1,\ldots,\boldsymbol{\alpha}_K)\)

Classiquement, on reparamétrise par

  • Pour chaque classe \(k\) et chaque variable \(j\) \[(\alpha_k^{\ j1},\ldots,\alpha_k^{\ jm_j}) \rightarrow (a_k^{\ j1},\ldots,a_k^{\ jm_j},\varepsilon_k^{\ j1},\ldots,\varepsilon_k^{\ jm_j})\] avec \[ a_k^{\ jh}=\left\{\begin{array}{l l}1 & \textrm{ si } h=\underset{h'=1,\ldots,m_j}{\mbox{argmax}}\ \alpha_k^{\ jh'}\\ 0 & \textrm{sinon}\end{array}\right. \] (\(h\)= modalité majoritaire pour variable \(j\) dans classe \(k\)) et

\[ \varepsilon_{k}^{\ jh}=\left\{\begin{array}{l l}1 - \alpha_k^{\ jh} & \textrm{ si } a_k^{\ jh}=1\\ \alpha_k^{\ jh} & \textrm{ si } a_k^{\ jh}=0\end{array}\right. \]

Par exemple \((0.3,0.6,0.1) \rightsquigarrow (0,1,0,\ \ 0.3,0.4,0.1)\)

La densité de mélange se réécrit alors \[f(.|\theta_K)=\underset{k=1}{\stackrel{K}{\sum}}\pi_k \underset{j=1}{\stackrel{p}{\prod}}\underset{h=1}{\stackrel{m_j}{\prod}} \left[\left(1-\varepsilon_k^{\,jh}\right)^{a_k^{\ jh}} \left(\varepsilon_k^{\,jh}\right)^{1-a_k^{\ jh}}\right]^{x_i^{jh}}\]

Selon les hypothèses faites sur les \(\varepsilon_{k}^{\ jh}\) et sur les proportions du mélange, on a 10 formes possibles (dans le même esprit que les 24 formes des mélanges gaussiens, cf cours).

Question : A l’aide de la fonction mixmodCluster de la librairie Rmixmod, déterminez une classification des données par modèles de mélange. Vous étudierez en particulier les probabilités conditionnelles d’appartenance.

library(Rmixmod)

resmixmod <- mixmodCluster(data = Data, nbCluster = 2:15, criterion = c("BIC", "ICL"),
    model = mixmodMultinomialModel("Binary_pk_Ekjh"))

# Graphe des critères de sélection
K <- NULL
BIC <- NULL
ICL <- NULL
for (k in 1:length(resmixmod@results)) {
    K <- c(K, resmixmod@results[[k]]@nbCluster)
    BIC <- c(BIC, resmixmod@results[[k]]@criterionValue[1])
    ICL <- c(ICL, resmixmod@results[[k]]@criterionValue[2])
}

## graphique à faire
...


# Etude de la classification retenue
df <- data.frame(proba = apply(resmixmod@bestResult@proba, 1, max), label = as.factor(apply(resmixmod@bestResult@proba,
    1, which.max)))
ggplot(df, aes(x = label, y = proba)) + geom_boxplot()
table(resmixmod@bestResult@partition)
## A completer

# Comparaison avec les autres classifications A completer

Correction :

library(Rmixmod)

resmixmod <- mixmodCluster(data = Data, nbCluster = 2:15, criterion = c("BIC", "ICL"),
    model = mixmodMultinomialModel("Binary_pk_Ekjh"))
# Graphe des critères de sélection
K <- NULL
BIC <- NULL
ICL <- NULL
for (k in 1:length(resmixmod@results)) {
    K <- c(K, resmixmod@results[[k]]@nbCluster)
    BIC <- c(BIC, resmixmod@results[[k]]@criterionValue[1])
    ICL <- c(ICL, resmixmod@results[[k]]@criterionValue[2])
}

aux <- sort.int(K, index.return = T)
df <- data.frame(K = K[aux$ix], BIC = BIC[aux$ix], ICL = ICL[aux$ix])
df <- melt(df, id = c("K"))
ggplot(df, aes(x = K, y = value, color = variable)) + geom_line()

# Etude de la classification retenue
df <- data.frame(proba = apply(resmixmod@bestResult@proba, 1, max), label = as.factor(apply(resmixmod@bestResult@proba,
    1, which.max)))
ggplot(df, aes(x = label, y = proba)) + geom_boxplot()

table(resmixmod@bestResult@partition)

  1   2   3   4   5   6   7 
 69 188 127 126 198  28  64 
barplotClus(resmixmod@bestResult@partition, Data, J = c(1:10))

fviz_mca_ind(resacm, habillage = as.factor(resmixmod@bestResult@partition), geom = "point") +
    labs(title = "Mixture")

# Comparaison avec Class
table(resmixmod@bestResult@partition, classesp)
   classesp
      e   p
  1  69   0
  2   0 188
  3   0 127
  4  77  49
  5 181  17
  6  24   4
  7  49  15
adjustedRandIndex(resmixmod@bestResult@partition, classesp)
[1] 0.2655892
# Comparaison avec Kmeans
table(resmixmod@bestResult@partition, kmeansclus[, Kaux])
   
      1   2   3   4   5   6   7   8
  1   0   0   0   0   0   0  69   0
  2   0 188   0   0   0   0   0   0
  3   0   0 127   0   0   0   0   0
  4   0   0   0   0   0  21 105   0
  5   0   0   0   0   0 195   3   0
  6  17   2   0   4   0   5   0   0
  7   0   0   0  41   6   0   0  17
adjustedRandIndex(resmixmod@bestResult@partition, kmeansclus[, Kaux])
[1] 0.8255953